In this sub-vignette we present the analysis and source code for figure 5. This sub-vignette can be built along with all other sub-vignettes, by running CLLCytokineScreen2021.Rmd.

1 Set up

Load Libraries

library(patchwork)
library(ggplot2)
library(data.table)
library(magrittr)
library(ggplot2)
library(pheatmap)
library(tidyr)
library(dplyr)

Set plot directory

plotDir = ifelse(exists(".standalone"), "", "../../inst/figs/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)

2 Load data

Raw

#df: tibble containing all screening viability data
load( "../../data/df.RData")
#from tsvs
df <- read.table(file= "../../inst/extdata/df.txt",header = TRUE) %>% as_tibble()

df$Cytokine <- factor(df$Cytokine, levels= c("No Cytokine",
                                              "Resiquimod", 
                                              "IL-4", 
                                              "TGF-b1",
                                              "IL-1b",
                                              "Interferon gamma",
                                              "SDF-1a",
                                              "sCD40L" ,
                                              "sCD40L+IL-4",
                                              "soluble anti-IgM", 
                                              "CpG ODN",
                                              "IL-6",
                                              "IL-10",
                                              "IL-21",
                                              "HS-5 CM", 
                                              "IL-15",
                                              "BAFF",
                                              "IL-2" ))

df$Drug <- factor(df$Drug, levels =c("DMSO",
                                     "BAY-11-7085",
                                     "Everolimus",
                                     "Fludarabine",
                                     "I-BET 762",
                                     "Ibrutinib","Idelalisib",
                                     "Luminespib",
                                     "Nutlin-3a",
                                     "PRT062607",
                                     "Pyridone-6",
                                     "Ralimetinib",
                                     "Selumetinib"))

Process data

# List of Cytokines, Drugs and treatments 
  ##Drugs
  thedrugs <- unique(df$Drug) 

  ##Cytokines
  thecytokines <- unique(df$Cytokine) 

  #Treatments
  #drugs
  drugtreatments <- list()
  drugtreatments <- 
    lapply(thedrugs, function(x){
      paste("treatment_drug", x, sep='')
      })
  
  
  #cytokines
  cytokinetreatments <- list()
  cytokinetreatments <- 
    lapply(thecytokines, function(y){
      paste("treatment_cytokine", y, sep='')
      })


# Tidy up data frame to use in linear model
df <- 
  dplyr::filter(df,
                #use high concentrations only
                Drug_Concentration %in% c("High","None")) %>% 
  #select columns for linear model
  dplyr::select(PatientID, Drug, Cytokine, Log) %>%  
        
  #rename columns 
  plyr::rename(c("Drug"="treatment_drug", "Cytokine"="treatment_cytokine", "Log"="Viability"))

3 Set aesthetics

source("../../R/themes_colors.R")

4 Pre - work

Fit the linear model

# define the base level per treatment_type as the "no"-treatment
df$treatment_drug <- as.factor(df$treatment_drug)
df$treatment_cytokine <- as.factor(df$treatment_cytokine)

df$treatment_drug %<>% relevel("DMSO")
df$treatment_cytokine %<>% relevel("No Cytokine")


# inspect design matrix of interaction model
# model.matrix(~treatment_drug + treatment_cytokine + treatment_drug:treatment_cytokine, df)

# fit model with interaction
fit <- lm(Viability ~ treatment_drug * treatment_cytokine, df)

Extract p values

#extract p values

pvalues <- summary(fit)$coefficients[,4]

#order as a dataframe
pvaldf <- as.data.frame(as.matrix(pvalues)) %>% 
          setDT(keep.rownames = TRUE)
          
#rename columns
colnames(pvaldf) <- c("treatment", "pvalue")

#filter out single agent treatments
singletreatments <- c(drugtreatments, cytokinetreatments, "(Intercept)")

pvaldf <- dplyr::filter(pvaldf, !treatment %in% singletreatments)


pvaldf <- 
  pvaldf %>% 
  #remove "treatment" prefix
  mutate_at(vars(treatment), funs(as.character(gsub("treatment_drug", "", .)))) %>% 
  mutate_at(vars(treatment), funs(as.character(gsub("treatment_cytokine", "", .))))
## Warning: `funs()` was deprecated in dplyr 0.8.0.
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
  #split drug:cytokine into two columns
pvaldf <-   
  data.frame(pvaldf, do.call(rbind, strsplit(pvaldf$treatment, split = ":", fixed = TRUE)))

#renove treatment column
pvaldf <- pvaldf[,c("pvalue","X1", "X2")]

#rename columns
colnames(pvaldf) <- c("pvalue","drug", "Cytokine")

5 Plot Figures

5.1 Fig 5A

Dummy plots

dummydata <-
  as.data.frame(list(c(1,2,3,4),
                     c(0,0,0,0),
                     c(-0.6,-0.6,0,-0.1), 
                     c(0.6,0.6,0.1,-0.1), 
                     c(0.5,-0.5,0.8,-0.8)), 
                col.names=c("plot",
                            "coord1", 
                            "coord2", 
                            "coord3", 
                            "coord4")) %>% 
  as_tibble()%>% 
  
  mutate(expected = coord2 + coord3, measured = coord4) %>% 
  
  pivot_longer(cols = coord1:measured, names_to = c("plot"), values_to="y.value", names_repair="unique" ) %>% 
  
  dplyr::rename(Plot = plot...1, x.value = plot...2) %>% 
  
  mutate(x.value = factor(x.value, levels = c("coord1", "coord2", "coord3", "coord4", "expected", "measured"))) %>% 
  
  #add annotation of interaction type
  mutate(Plot = factor(Plot, 
                       labels = c("stimulus inhibits drug", 
                                  "drug inhibits stimulus", 
                                  "synergistic pro-survival effect", 
                                  "synergistic toxicity")))
## New names:
## * plot -> plot...1
## * plot -> plot...2
gg_dummy =
  #for each interaction type, make a dummy plot
  lapply(unique(dummydata$Plot), function(i){
    
    tmp <- dplyr::filter(dummydata, Plot==i )
    
    #get value for interaction
    interaction = 
      unlist(dplyr::select(dplyr::filter(tmp, x.value=="measured"), y.value)) - 
      unlist(dplyr::select(dplyr::filter(tmp,x.value=="expected"), y.value))
    
    #plot dummy plot
    gg_dummy <- 
      ggplot(dplyr::filter(tmp, grepl('coord', x.value)), 
             aes(x=x.value, y=y.value, group=Plot))+
  
  geom_hline(yintercept = 0, linetype=2, color="black")+
  geom_point(colour=ifelse(interaction > 0, "#A6093D","#003DA5"), size=1)+
  geom_line(colour=ifelse(interaction > 0, "#A6093D","#003DA5"), size=1)+
  
  geom_segment(data=dplyr::filter(tmp, grepl('expected', x.value)), aes(x = 3.5, y = y.value, xend = 4.5, yend = y.value), color="blue", size=2) +
  geom_segment(data=dplyr::filter(tmp, grepl('measured', x.value)), aes(x = 3.5, y = y.value, xend = 4.5, yend = y.value), color="black", size=2)+
  
  theme_bw() +
  t1 +
  xlab("") + 
  ylab("") +
  coord_cartesian(ylim = c(-1, 1)) +
  scale_x_discrete(labels=c( "coord1"="DMSO", "coord2"="Drug","coord3"="Stimulus","coord4"="Drug +\nStimulus")) +
    theme(axis.ticks.x= element_blank(),axis.ticks.y = element_blank(),
          plot.tag=element_text(size = 20)) +
      
  scale_y_continuous(breaks = c(-1,0,1), sec.axis = 
                       sec_axis(trans=~exp(1)^.*100, 
                                breaks=c(50,100,200), 
                                labels=c("50%", "100%", "200%")))

})


#construct plot

Fig5A <- wrap_elements(
  
  gg_dummy[[1]] + 
  labs(title=expression(paste("Positive ", beta["int"] ))) +
  ylab("Antagonism")+

  gg_dummy[[2]] +
  labs(title=expression(paste("Negative ", beta["int"] )))+

  gg_dummy[[3]] +
  ylab("Synergism") +
    
  gg_dummy[[4]]+
  
  plot_layout(nrow=2) +
  plot_annotation(tag_levels = "I")
)

Fig5A

5.2 Fig 5B

Barplot quantifying interaction types

fit_as_df <- 
  fit$coefficients %>% 
  as.matrix() %>% 
  as.data.frame() %>%
  setDT(keep.rownames = TRUE) %>% 
  as_tibble() %>% 
  dplyr::rename(Condition = rn, comb_coefficient = V1) %>% 
  dplyr::filter(grepl("drug", Condition)&grepl("cytokine", Condition)) %>% 
  mutate(Condition = gsub("treatment_drug", "", Condition)) %>% 
  mutate(Condition = gsub("treatment_cytokine", "", Condition)) %>% 
  separate(Condition, c("Drug", "Cytokine"), sep=":")

drug_fit_as_df <- 
  fit$coefficients %>% 
  as.matrix() %>%
  as.data.frame() %>%
  setDT(keep.rownames = TRUE) %>% 
  as_tibble() %>% 
  dplyr::rename(Drug = rn, drug_coefficent=V1) %>% 
  dplyr::filter(grepl("drug", Drug)&!grepl("cytokine", Drug)) %>% 
  mutate(Drug = gsub("treatment_drug", "", Drug))


cyt_fit_as_df <-
  fit$coefficients %>% 
  as.matrix() %>% 
  as.data.frame() %>%
  setDT(keep.rownames = TRUE) %>% 
  as_tibble() %>% 
  dplyr::rename(Cytokine=rn, cyt_coefficent=V1) %>% 
  dplyr::filter(!grepl("drug", Cytokine)&grepl("cytokine", Cytokine)) %>% 
  mutate(Cytokine = gsub("treatment_cytokine", "", Cytokine))


fit_combinations <-
  left_join(fit_as_df, drug_fit_as_df, by = "Drug") %>% 
  left_join(cyt_fit_as_df, by = "Cytokine") %>% 
  mutate(predicted_coefficient=cyt_coefficent+drug_coefficent) %>% 
  mutate(obs_coef=predicted_coefficient+comb_coefficient) %>%
  left_join(mutate(pvaldf,drug=as.character(drug),Cytokine=as.character(Cytokine)), by=c("Drug"="drug", "Cytokine"="Cytokine")) %>% 
  dplyr::filter(pvalue<=0.05) %>%
  mutate(Category_Symbol=case_when(
    comb_coefficient>=0&(obs_coef<cyt_coefficent|obs_coef<drug_coefficent)~ "I",  
    comb_coefficient>=0&(obs_coef>cyt_coefficent&obs_coef>drug_coefficent)~ "III",
    comb_coefficient<0&(obs_coef>cyt_coefficent|obs_coef>drug_coefficent)~ "II",  
    comb_coefficient<0&(obs_coef<cyt_coefficent&obs_coef<drug_coefficent)~ "IV"
    )) %>% 
  mutate(Category=case_when(
    Category_Symbol=="I"~ "Positive Antagonistic",  
    Category_Symbol=="III"~ "Positive Synergistic",
    Category_Symbol=="II"~ "Negative Antagonistic",  
    Category_Symbol=="IV"~ "Negative Synergistic"
    )) %>% 
  mutate( Category=factor(Category, levels=c("Positive Antagonistic","Negative Antagonistic","Positive Synergistic","Negative Synergistic")))
  

Fig5B <-
  fit_combinations %>% 
  dplyr::group_by(Category) %>% 
  tally() %>% 
  
ggplot(aes(x=Category, y=n, fill=Category))+
  geom_bar(width = 1, stat = "identity")+
  t1+
  guides(fill="none")+
  xlab("")+ ylab("")+
  scale_fill_manual(values = c("#A6093D","#003DA5","#A6093D","#003DA5"))

Fig5B

5.3 Fig 5C

Heatmap of significant drug : stimulus interaction coefficients
Get interaction coefficients matrix for plotting heatmap

#extract beta interaction values
betavalues <- summary(fit)$coefficients[,1]

#create matrix of beta values 
betadf <- as.data.frame(as.matrix(betavalues)) %>%
          setDT(keep.rownames = TRUE)

#rename columns
colnames(betadf) <- c("treatment", "betavalue")

#remove beta values for drugs and cytokines alone
betadf <- dplyr::filter(betadf, !treatment %in% singletreatments)

#remove treatment prefix
betadf <- betadf %>% 
  mutate_at(vars(treatment), funs(as.character(gsub("treatment_drug", "", .)))) %>% 
  mutate_at(vars(treatment), funs(as.character(gsub("treatment_cytokine", "", .))))

#split drug:cytokine by colon
betadf <- data.frame(betadf, do.call(rbind, strsplit(betadf$treatment, split = ":", fixed = TRUE)))

#select columns of interest
betadf <- betadf[,c("betavalue","X1", "X2")]

#rename columns 
colnames(betadf) <- c("betavalue","drug", "Cytokine")
#set significance
a <- 0.05

#make a matrix of beta values, where  beta = 0 if corresponding p val is > than significance threshold 
pvalmat <- xtabs(pvalue~drug+Cytokine, data=pvaldf)
bvalmat <- xtabs(betavalue~drug+Cytokine, data=betadf)

#check in same order
bvalmat[,colnames(pvalmat)] 
##              Cytokine
## drug                   BAFF       CpG ODN       HS-5 CM         IL-10
##   BAY-11-7085 -5.604631e-02 -3.181140e-01  1.492971e-03 -1.575975e-02
##   Everolimus   3.655902e-02  1.497761e-02  5.377161e-02  7.948169e-02
##   Fludarabine -2.337342e-02  1.903315e-01  6.711677e-02 -3.913832e-02
##   I-BET 762    4.216039e-02 -2.331726e-01 -7.116831e-02 -2.652204e-02
##   Ibrutinib   -2.625863e-02 -2.269246e-01  7.477386e-02  2.731021e-02
##   Idelalisib   8.275484e-03 -2.652441e-01  3.386327e-02  3.058002e-02
##   Luminespib  -2.482555e-01 -2.591311e-01 -3.940346e-01 -6.458359e-02
##   Nutlin-3a   -1.258362e-03  1.577514e-01  7.700330e-02  9.299536e-02
##   PRT062607   -6.653378e-02 -8.644230e-02  1.557073e-02 -8.013149e-03
##   Pyridone-6  -8.023060e-04 -8.100430e-03  2.925389e-02  3.243689e-02
##   Ralimetinib -3.093962e-02 -1.717127e-01  1.858811e-02 -1.073542e-02
##   Selumetinib  5.947874e-03 -1.304805e-01  7.060345e-02  2.145764e-02
##              Cytokine
## drug                  IL-15         IL-1b          IL-2         IL-21
##   BAY-11-7085 -5.816512e-03  2.152596e-02 -2.134484e-02  2.072687e-03
##   Everolimus   2.725033e-02  3.500697e-02 -2.677401e-03  3.495098e-02
##   Fludarabine -4.310581e-02  3.505833e-02 -2.214672e-02 -4.860353e-02
##   I-BET 762    9.141791e-04  2.740045e-03 -2.907140e-02 -2.523721e-02
##   Ibrutinib    5.752202e-02  1.849707e-03  4.678019e-02  3.710308e-02
##   Idelalisib   5.658295e-02  1.936971e-02  4.584077e-02  3.477902e-02
##   Luminespib  -6.061613e-02 -1.263584e-01 -6.993392e-02 -7.742445e-02
##   Nutlin-3a    1.046024e-01  2.482328e-02  1.463441e-01  7.912615e-02
##   PRT062607   -3.957221e-03 -2.216149e-02  5.916628e-03  8.482059e-03
##   Pyridone-6   5.308392e-02  1.514414e-02  3.040022e-02 -7.237789e-04
##   Ralimetinib  1.202065e-02 -2.568524e-03 -1.468196e-02 -4.467606e-02
##   Selumetinib  5.942633e-02  1.780087e-02  3.092410e-02  2.662341e-02
##              Cytokine
## drug                   IL-4          IL-6 Interferon gamma    Resiquimod
##   BAY-11-7085  2.235761e-02  4.554573e-02     5.434015e-02 -3.639915e-01
##   Everolimus   1.332949e-01  6.817609e-03     5.469766e-02  2.765104e-02
##   Fludarabine  3.953363e-01 -2.469675e-02     4.886592e-02  2.265297e-01
##   I-BET 762   -2.665207e-02 -2.262955e-02     1.357931e-02 -3.265356e-01
##   Ibrutinib    1.836975e-01  6.642272e-02     1.101352e-01 -2.509692e-01
##   Idelalisib   1.657969e-01  6.363235e-02     1.171099e-01 -2.064261e-01
##   Luminespib  -4.429811e-02 -3.220023e-03    -2.221338e-02 -7.116024e-02
##   Nutlin-3a    5.337703e-01  1.136707e-01     1.433726e-01  1.426722e-01
##   PRT062607    1.388865e-01 -1.097500e-02     2.337567e-02 -3.612143e-01
##   Pyridone-6  -2.305619e-02  6.847247e-02     6.271102e-02  1.349796e-03
##   Ralimetinib  7.075108e-02  2.282944e-03     2.190150e-01 -1.118603e-01
##   Selumetinib  1.174186e-01  6.322721e-02     8.349273e-02 -1.148979e-01
##              Cytokine
## drug                 sCD40L   sCD40L+IL-4        SDF-1a soluble anti-IgM
##   BAY-11-7085  2.049505e-02 -1.619583e-01  2.831627e-02    -6.597054e-03
##   Everolimus   5.178276e-02  1.386887e-01 -2.677009e-03     5.541420e-04
##   Fludarabine  1.020548e-02  3.687631e-01  5.075778e-02     3.435734e-02
##   I-BET 762    4.068479e-02  6.576262e-03 -4.944505e-03    -2.502231e-02
##   Ibrutinib   -4.740237e-03  8.873182e-02  2.038604e-02     1.260542e-02
##   Idelalisib   2.027841e-02  8.430112e-02 -2.252409e-05     2.338468e-02
##   Luminespib  -1.543305e-01 -3.482382e-01 -8.004006e-02    -4.275447e-01
##   Nutlin-3a    1.834396e-02  4.643749e-01  1.258999e-01     5.912671e-03
##   PRT062607   -8.838768e-02  1.457566e-01 -2.877775e-02    -6.326385e-02
##   Pyridone-6   3.934846e-02 -1.855970e-01  3.145554e-02    -6.082148e-03
##   Ralimetinib -1.691833e-04 -3.878371e-02  1.575747e-03    -2.273022e-02
##   Selumetinib  5.122307e-02  4.874905e-02 -2.510831e-02     2.421239e-02
##              Cytokine
## drug                 TGF-b1
##   BAY-11-7085 -4.509014e-03
##   Everolimus  -8.310092e-03
##   Fludarabine -4.636054e-02
##   I-BET 762   -9.396045e-02
##   Ibrutinib    1.958423e-02
##   Idelalisib   9.050594e-03
##   Luminespib  -8.478118e-02
##   Nutlin-3a   -1.529464e-02
##   PRT062607    2.369375e-03
##   Pyridone-6   3.081690e-02
##   Ralimetinib  1.357279e-02
##   Selumetinib  2.234705e-02
bvalmat[rownames(pvalmat),]
##              Cytokine
## drug                   BAFF       CpG ODN       HS-5 CM         IL-10
##   BAY-11-7085 -5.604631e-02 -3.181140e-01  1.492971e-03 -1.575975e-02
##   Everolimus   3.655902e-02  1.497761e-02  5.377161e-02  7.948169e-02
##   Fludarabine -2.337342e-02  1.903315e-01  6.711677e-02 -3.913832e-02
##   I-BET 762    4.216039e-02 -2.331726e-01 -7.116831e-02 -2.652204e-02
##   Ibrutinib   -2.625863e-02 -2.269246e-01  7.477386e-02  2.731021e-02
##   Idelalisib   8.275484e-03 -2.652441e-01  3.386327e-02  3.058002e-02
##   Luminespib  -2.482555e-01 -2.591311e-01 -3.940346e-01 -6.458359e-02
##   Nutlin-3a   -1.258362e-03  1.577514e-01  7.700330e-02  9.299536e-02
##   PRT062607   -6.653378e-02 -8.644230e-02  1.557073e-02 -8.013149e-03
##   Pyridone-6  -8.023060e-04 -8.100430e-03  2.925389e-02  3.243689e-02
##   Ralimetinib -3.093962e-02 -1.717127e-01  1.858811e-02 -1.073542e-02
##   Selumetinib  5.947874e-03 -1.304805e-01  7.060345e-02  2.145764e-02
##              Cytokine
## drug                  IL-15         IL-1b          IL-2         IL-21
##   BAY-11-7085 -5.816512e-03  2.152596e-02 -2.134484e-02  2.072687e-03
##   Everolimus   2.725033e-02  3.500697e-02 -2.677401e-03  3.495098e-02
##   Fludarabine -4.310581e-02  3.505833e-02 -2.214672e-02 -4.860353e-02
##   I-BET 762    9.141791e-04  2.740045e-03 -2.907140e-02 -2.523721e-02
##   Ibrutinib    5.752202e-02  1.849707e-03  4.678019e-02  3.710308e-02
##   Idelalisib   5.658295e-02  1.936971e-02  4.584077e-02  3.477902e-02
##   Luminespib  -6.061613e-02 -1.263584e-01 -6.993392e-02 -7.742445e-02
##   Nutlin-3a    1.046024e-01  2.482328e-02  1.463441e-01  7.912615e-02
##   PRT062607   -3.957221e-03 -2.216149e-02  5.916628e-03  8.482059e-03
##   Pyridone-6   5.308392e-02  1.514414e-02  3.040022e-02 -7.237789e-04
##   Ralimetinib  1.202065e-02 -2.568524e-03 -1.468196e-02 -4.467606e-02
##   Selumetinib  5.942633e-02  1.780087e-02  3.092410e-02  2.662341e-02
##              Cytokine
## drug                   IL-4          IL-6 Interferon gamma    Resiquimod
##   BAY-11-7085  2.235761e-02  4.554573e-02     5.434015e-02 -3.639915e-01
##   Everolimus   1.332949e-01  6.817609e-03     5.469766e-02  2.765104e-02
##   Fludarabine  3.953363e-01 -2.469675e-02     4.886592e-02  2.265297e-01
##   I-BET 762   -2.665207e-02 -2.262955e-02     1.357931e-02 -3.265356e-01
##   Ibrutinib    1.836975e-01  6.642272e-02     1.101352e-01 -2.509692e-01
##   Idelalisib   1.657969e-01  6.363235e-02     1.171099e-01 -2.064261e-01
##   Luminespib  -4.429811e-02 -3.220023e-03    -2.221338e-02 -7.116024e-02
##   Nutlin-3a    5.337703e-01  1.136707e-01     1.433726e-01  1.426722e-01
##   PRT062607    1.388865e-01 -1.097500e-02     2.337567e-02 -3.612143e-01
##   Pyridone-6  -2.305619e-02  6.847247e-02     6.271102e-02  1.349796e-03
##   Ralimetinib  7.075108e-02  2.282944e-03     2.190150e-01 -1.118603e-01
##   Selumetinib  1.174186e-01  6.322721e-02     8.349273e-02 -1.148979e-01
##              Cytokine
## drug                 sCD40L   sCD40L+IL-4        SDF-1a soluble anti-IgM
##   BAY-11-7085  2.049505e-02 -1.619583e-01  2.831627e-02    -6.597054e-03
##   Everolimus   5.178276e-02  1.386887e-01 -2.677009e-03     5.541420e-04
##   Fludarabine  1.020548e-02  3.687631e-01  5.075778e-02     3.435734e-02
##   I-BET 762    4.068479e-02  6.576262e-03 -4.944505e-03    -2.502231e-02
##   Ibrutinib   -4.740237e-03  8.873182e-02  2.038604e-02     1.260542e-02
##   Idelalisib   2.027841e-02  8.430112e-02 -2.252409e-05     2.338468e-02
##   Luminespib  -1.543305e-01 -3.482382e-01 -8.004006e-02    -4.275447e-01
##   Nutlin-3a    1.834396e-02  4.643749e-01  1.258999e-01     5.912671e-03
##   PRT062607   -8.838768e-02  1.457566e-01 -2.877775e-02    -6.326385e-02
##   Pyridone-6   3.934846e-02 -1.855970e-01  3.145554e-02    -6.082148e-03
##   Ralimetinib -1.691833e-04 -3.878371e-02  1.575747e-03    -2.273022e-02
##   Selumetinib  5.122307e-02  4.874905e-02 -2.510831e-02     2.421239e-02
##              Cytokine
## drug                 TGF-b1
##   BAY-11-7085 -4.509014e-03
##   Everolimus  -8.310092e-03
##   Fludarabine -4.636054e-02
##   I-BET 762   -9.396045e-02
##   Ibrutinib    1.958423e-02
##   Idelalisib   9.050594e-03
##   Luminespib  -8.478118e-02
##   Nutlin-3a   -1.529464e-02
##   PRT062607    2.369375e-03
##   Pyridone-6   3.081690e-02
##   Ralimetinib  1.357279e-02
##   Selumetinib  2.234705e-02
bvalmat[pvalmat >= a] = 0

#transform matrix 
bvalmat <- t(bvalmat)

tree <- 
  pheatmap(bvalmat, silent = TRUE)


cyt_order <- rownames(bvalmat[tree$tree_row[["order"]],])

drug_order <- colnames(bvalmat[,tree$tree_col[["order"]]])


Fig5C <-
  #append p values
  left_join(betadf, pvaldf, by = c("drug", "Cytokine")) %>% 
  #make beta value 0 if p value is < sig
  mutate(betavalue = ifelse(pvalue<a, betavalue, 0)) %>% 
  #put drugs in order of clustered heatmap
  mutate(drug=factor(drug, levels = drug_order)) %>% 
  
  #put stimuli in order of clustered heatmap
  mutate(Cytokine=factor(Cytokine, levels = rev(cyt_order))) %>% 
  
  #make heatmap with ggplot
  ggplot(aes(x=drug, y=Cytokine))+
  
  geom_tile(aes(fill=betavalue),color = "grey")+
  #set colour scale
  scale_fill_gradientn(colors=c("#003DA5",  "white",  "#A6093D"), limits=c(-0.6,0.6))+
  #set theme
  t2+
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        axis.ticks.x =element_blank(),
        panel.background = element_blank(), 
        axis.line.x = element_blank(),
        axis.line.y = element_blank(),
        legend.key.height=unit(2.5, "cm"))+
  labs(fill = expression(paste(beta["int"], "-value")))+
  #add category
  geom_text(data = fit_combinations, aes(x=Drug, label=Category_Symbol) )+
  scale_y_discrete(labels=c("TGF-b1"="TGF-\u03B21", "sCD40L+IL-4"="sCD40L + IL4", "IL-1b"="IL1\u03B2", "IL-4"="IL4", "IL-6"="IL6","IL-15"="IL15","IL-10"="IL10", "IL-21"="IL21","IL-2"="IL2", "Interferon gamma"= "Interferon \u03B3", "SDF-1a"="SDF-1\u03B1"))
  
Fig5C

5.4 Fig 5D - I

Preparation: Line plots of examples of drug and stimulus combinations that show an interaction

#get lists of drugs and cytokines without baseline treatments
thedrugs.only <- thedrugs %>% setdiff("DMSO") %>% sort()
thecytokines.only <- thecytokines %>% setdiff("No Cytokine")%>% sort()

# get a list of interactions for which p value of beta interaction is significant
sign_conditions <-
  pvaldf %>% 
  dplyr::filter(pvalue<=0.05) %>% 
  mutate(condition = paste0("treatment_drug:", 
                            drug, 
                            ":treatment_cytokine:",
                            Cytokine))


gg <- vector(mode="list", length=length(sign_conditions$condition))
names(gg) <- sign_conditions$condition



plotList <- 
lapply(names(gg), function(i){
  
  #preparations for the interaction plot
  drug <- dplyr::filter(sign_conditions, condition ==i)$drug
  drugTreat <- paste("treatment_drug", drug, sep='')
      
  cyt <- dplyr::filter(sign_conditions, condition ==i)$Cytokine
  cytokineTreat <-paste("treatment_cytokine", cyt, sep='')
      
      
  Category <- fit_combinations[which(fit_combinations$Drug==drug &
                               fit_combinations$Cytokine==cyt),]$Category
  
  #extract coefficients to plot

  ## baseline effect
  baseline <- fit$coefficients["(Intercept)"]

  ## single treatment effect for type_1 = "a"
  single_1 <- fit$coefficients[drugTreat]

    ## single treatment effect for type_2 = "b"
    single_2 <- fit$coefficients[cytokineTreat]

    # interaction Viability
    interaction <- fit$coefficients[paste(drugTreat, cytokineTreat, sep=':')]

      
    ## observed single treatment values
    observed_single_1 <- baseline + single_1 
    observed_single_2 <- baseline + single_2

    ## observed double treatment value
    observed_double <- baseline + single_1 + single_2 + interaction

    ## predicted double treatment value (based only on single effects)
    predicted_double <- baseline + single_1 + single_2

    #prepare labels and aesthetics for plotting
    # y_label <- case_when(cyt=="Interferon gamma"~"Interferon\ngamma",
    #                      cyt=="soluble anti-IgM"~"soluble\nanti-IgM",
    #                      TRUE~cyt)

    xlabels <- c("DMSO", drug, 
                 case_when(cyt=="IL-4"~"IL4",
                           cyt=="Interferon gamma"~"Interferon \u03B3",
                           cyt=="sCD40L+IL-4"~"sCD40L + IL4",
                           cyt=="soluble anti-IgM"~"soluble\nanti-IgM",
                           TRUE~cyt),
                paste(drug, " +\n",
                       case_when(cyt=="IL-4"~"IL4",
                           cyt=="Interferon gamma"~"Interferon \u03B3",
                           cyt=="sCD40L+IL-4"~"sCD40L + IL4",
                           TRUE~cyt),
                sep=""))
      
    segment_size <- 2
    
    
    plotTab <-
      df %>%
      
      dplyr::filter(treatment_drug %in% c("DMSO" , drug)) %>%

      dplyr::filter(treatment_cytokine %in% c("No Cytokine", cyt)) %>%

      dplyr::mutate(treatment_combination = 
                  #if baseline treatment
                  ifelse(treatment_drug=="DMSO" & 
                         treatment_cytokine=="No Cytokine",
                         "DMSO",
                         #if drug only treatment
                         ifelse(treatment_drug == drug & 
                                treatment_cytokine =="No Cytokine",
                                paste0("Drug=", drug),
                                #if cytokine only treatment
                                ifelse(treatment_drug == "DMSO" & 
                                       treatment_cytokine == cyt,
                                       paste0("Cytokine=", cyt),
                                       #otherwise combinatorial treatment
                                       paste0("Drug=", drug, "\nCytokine=", cyt))))) %>%
  #make treatment_combination as a factor
  mutate(treatment_combination = factor(treatment_combination, 
                                        levels=c("DMSO", 
                                                 paste0("Drug=", drug),
                                                 paste0("Cytokine=", cyt),
                                                 paste0("Drug=", drug, "\nCytokine=", cyt)
                                                 )))
  #make plot
        
        

  ggplot(plotTab, aes(treatment_combination, Viability, group=PatientID)) +
  
  geom_hline(yintercept = 0, linetype=2, color="black")+

  lemon::geom_pointline( size=1,   
                         colour=ifelse(interaction > 0,
                                       "#A6093D","#003DA5"), alpha=0.3) +
  #add predicted viability with control treatment
  geom_segment( size=segment_size,
                aes(x=1-.2, xend=1+.2, y=baseline, yend=baseline),
                colour = "black") +
  #add predicted viability with drug only treatment
  geom_segment( size=segment_size,
                aes(x=2-.2, xend=2+.2, y=observed_single_1,yend=observed_single_1),
                colour="black") +
  #add predicted viability with cytokine only treatment
  geom_segment( size=segment_size,
                aes(x=3-.2, xend=3+.2, y=observed_single_2, yend=observed_single_2),
                colour="black") +
  #add predicted viability with both treatments
  geom_segment( size=segment_size,
                aes(x=4-.2, xend=4+.2, y=observed_double, yend=observed_double),
                colour="black") +
  #add predicted viability with both treatments, without interaction
  geom_segment( size=segment_size,
                aes(x=4-.2, xend=4+.2, y=predicted_double, 
                    yend=predicted_double), color="blue") + 
 
   #add category and combination as title
    
    
    
  ggtitle(paste( Category,": ", drug, " & ", 
                 case_when(cyt=="IL-4"~"IL4",
                           cyt=="Interferon gamma"~"Interferon \u03B3",
                           cyt=="sCD40L+IL-4"~"sCD40L + IL4",
                           TRUE~cyt), sep="")) + 
  
  #add theme
  t2 + theme(axis.text.x = element_text(size=fontsize+5)) +
    theme(plot.tag=element_text(size = 30))+
  scale_x_discrete(labels = xlabels) +
  ylab("Log (Viability)") + 
  xlab("")


  } ) 

names(plotList) <- names(gg)

5.5 Fig 5D

wrap_plots(plotList$`treatment_drug:Ibrutinib:treatment_cytokine:IL-4`)

5.6 Fig 5E

wrap_plots(plotList$`treatment_drug:Ibrutinib:treatment_cytokine:Interferon gamma`)

5.7 Fig 5F

wrap_plots(plotList$`treatment_drug:Pyridone-6:treatment_cytokine:sCD40L+IL-4`)

5.8 Fig 5G

wrap_plots(plotList$`treatment_drug:Ralimetinib:treatment_cytokine:Interferon gamma`)

5.9 Fig 5H

wrap_plots(plotList$`treatment_drug:Ibrutinib:treatment_cytokine:CpG ODN`)

5.10 Fig 5G

wrap_plots(plotList$`treatment_drug:Luminespib:treatment_cytokine:soluble anti-IgM`)

6 Arrange plots

design1 <-"
  1223
  4455
  6677
  8899
  "

Fig5A + theme(plot.tag=element_text(size = 30)) +

wrap_elements(Fig5B) + theme(plot.tag=element_text(size = 30)) +
  
wrap_elements(Fig5C) + theme(plot.tag=element_text(size = 30)) +
  
wrap_plots(plotList$`treatment_drug:Ibrutinib:treatment_cytokine:IL-4`)+
  
wrap_plots(plotList$`treatment_drug:Ibrutinib:treatment_cytokine:Interferon gamma`)+
  
wrap_plots(plotList$`treatment_drug:Pyridone-6:treatment_cytokine:sCD40L+IL-4`)+
  
wrap_plots(plotList$`treatment_drug:Ralimetinib:treatment_cytokine:Interferon gamma`)+
  
wrap_plots(plotList$`treatment_drug:Ibrutinib:treatment_cytokine:CpG ODN`)+
  
wrap_plots(plotList$`treatment_drug:Luminespib:treatment_cytokine:soluble anti-IgM`)+
  
plot_layout(design=design1,
            heights = c(1,0.5,0.5,0.5),
            widths = c(0.92,0.08,0.05,0.95)) +

  plot_annotation( tag_levels = "A", theme=theme(plot.title = element_text(face="plain", size=28, hjust = 0.5)))

7 Appendix

Sys.info()
##                                                                                             sysname 
##                                                                                            "Darwin" 
##                                                                                             release 
##                                                                                            "19.6.0" 
##                                                                                             version 
## "Darwin Kernel Version 19.6.0: Thu May  6 00:48:39 PDT 2021; root:xnu-6153.141.33~1/RELEASE_X86_64" 
##                                                                                            nodename 
##                                               "mac-huber20.academic.christs.cam.ac.uk.pool.embl.de" 
##                                                                                             machine 
##                                                                                            "x86_64" 
##                                                                                               login 
##                                                                                             "holly" 
##                                                                                                user 
##                                                                                             "holly" 
##                                                                                      effective_user 
##                                                                                             "holly"
sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] dplyr_1.0.7       tidyr_1.1.3       pheatmap_1.0.12   magrittr_2.0.1   
## [5] data.table_1.14.0 ggplot2_3.3.5     patchwork_1.1.1   BiocStyle_2.20.2 
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.6          highr_0.9           plyr_1.8.6         
##  [4] RColorBrewer_1.1-2  pillar_1.6.1        compiler_4.1.0     
##  [7] BiocManager_1.30.16 tools_4.1.0         digest_0.6.27      
## [10] lattice_0.20-44     evaluate_0.14       lifecycle_1.0.0    
## [13] tibble_3.1.2        gtable_0.3.0        pkgconfig_2.0.3    
## [16] rlang_0.4.11        DBI_1.1.1           magick_2.7.2       
## [19] yaml_2.2.1          lemon_0.4.5         xfun_0.24          
## [22] gridExtra_2.3       withr_2.4.2         stringr_1.4.0      
## [25] knitr_1.33          generics_0.1.0      vctrs_0.3.8        
## [28] grid_4.1.0          tidyselect_1.1.1    glue_1.4.2         
## [31] R6_2.5.0            fansi_0.5.0         rmarkdown_2.9      
## [34] bookdown_0.22.3     farver_2.1.0        purrr_0.3.4        
## [37] scales_1.1.1        ellipsis_0.3.2      htmltools_0.5.1.1  
## [40] assertthat_0.2.1    colorspace_2.0-2    labeling_0.4.2     
## [43] utf8_1.2.1          stringi_1.6.2       munsell_0.5.0      
## [46] crayon_1.4.1